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The scalar nonlinear Schrodinger (NLS) equation and a suitable discretization are well known in- 
tegrable systems which exhibit the phenomena of "effective" chaos. Vector generalizations of both the 
0^ , continuous and discrete system are discussed. Some attention is directed upon the issue of the integra- 

bility of a discrete version of the vector NLS equation. 

1 Introduction 

O ' 

c/3 ■ Integrable systems are usually thought to demonstrate regular temporal and spatial evolution. Philosoph- 

ically, it seems to be pleasing to separate integrable motion from chaotic dynamics. The latter being the 
"opposite" of the former. However, upon more careful investigation, there is not such a simple distinction. 
Many well known integrable systems are actually "on the verge" of exhibiting chaos whereby any perturba- 
tion: physical, or numerical-e.g. even roundoff effects, can push them "over the edge" and thereby excite 
chaotic dynamics. In one dimension, the classical pendulum is such an example. And in one space - one time 
dimension the well known integrable systems: the sine-Gordon and nonlinear Schrodinger equations with 
periodic boundary values are examples of integrable systems where complex and irregular dynamics can be 
excited in certain regions of phase space ^, The common feature in such situations is an initial 

value prescribed in the proximity of special regimes of phase space, referred to as homoclinic manifolds. 
Homoclinic manifolds are locations in phase space in which highly sensitive solution structures exist. We 
have, in earlier papers, discussed the notion of "effective" chaos ^(whereby any perturbation can excite 
irregular and unstable motion-both of which are manifestations of chaotic dynamics. 

Chaotic behavior in perturbed integrable systems is significant for two reasons: (i) Numerical schemes 
can be considered as perturbations of the original system. If the perturbation of the numerical scheme 
introduces chaos, the results of simulations will not reflect the dynamics of the unperturbed system, (ii) In 
physical applications, perturbations -from effects neglected in the derivation of the integrable equation can 
induce chaotic behavior behavior in the system. Even if this effective chaos is only present in certain regions 
of phase space, it is, in general, impossible to determine a priori where these regions are. 

Finite-difference discretizations of integrable PDE's can be considered both as perturbations of of the 
original PDE and as the basis for numerical schemes. Here, we will consider spatial finite-difference schemes. 
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These finite-difference discretizations can be integrated in time by well-known, highly-accurate methods for 
integrating ordinary differential equations. Numerical schemes based on finite differences are useful because, 
unlike spectral schemes, finite- difference schemes do not require special boundary conditions. 
We begin by first discussing the scalar nonlinear Schrodinger equation (NLS): 



iq t = Qxx + 2q\q\ 2 . (f) 

NLS is integrable via the inverse scattering transform (1ST) and has significant applications in physics, such 
as the evolution of small amplitude slowly varying wave packets in: deep water, nonlinear optics and plasma 
physics (see e.g. [||). 

Many discretizations of NLS have physical applications as discrete systems (see e.g. ||, [T(], [l3| [HJ). 
In particular, consider two conservative (in fact, Hamiltonian) finite-difference discretizations of NLS, the 
diagonal discrete NLS (DDNLS), 

i^-q n = -r^(q n -i - 2q n + q n +i) + 2\q n \ 2 q n (2) 
dt hr 

and the integrable discrete NLS (IDNLS), 



i— q n = J^(ln-i - 2q n + q n +i) + \qn\ 2 {q n +\ + q n -i)- (3) 

These two systems differ only in the discretization of the nonlinear term, yet they have very different 
properties. While both these discretizations retain the Hamiltonian structure of the PDE, IDNLS (||) is 
integrable via the 1ST Q] while DDNLS is not. Differences in the dynamics of these two systems shed light 
on the role of integrability in the development of chaos in perturbations. In fact, these two schemes can 
exhibit very different responses to the same initial data. In the case of periodic boundary conditions, the 
onset and structure of "effective" chaos is very different between (g) and (||) while, for the infinite lattice 
with a rapidly-decaying boundary condition, (||) has solitons and (||) does not. The affect of non-conservative 
perturbations to NLS has also been considered in the literature (e. g. Jl?]]), but we will not do so here. 
The vector generalization of NLS, 



«q* = Qxx + 2 ||q|| 2 q (4) 

where q is an V-component vector and ||-|| denotes the vector norm, can be integrated via the 1ST and has 
soliton solutions |l(| (actually, in |ll| only the case N = 2 was studied in detail; however the extension to 
the 7V-th order vector system is straightforward). 

Vector NLS (eq. hereafter referred to as VNLS) is applicable, physically, under the similar conditions 
as NLS (§) with the generalization that the field can have more than one nontrivial component. For example, 
in optical fibers, the two components of the second-order (N = 2) VNLS correspond to components of the 
electric field transverse to the direction of wave propagation. These components of the transverse field 
compose a basis of the polarization states (we note that, in general, a non-integrable variation of VNLS is 
the appropriate model for optical fibers [Fst , however, in certain circumstances (Q) is indeed the appropriate 
equation (T^ 0). 

Although certain perturbations of the integrable VNLS have been investigated Jlj| |2l| , the question of 
effective chaos in VNLS has not been studied. Discretizations of VNLS are natural candidates for the study 
of chaos induced by perturbations. The research on effective chaos for NLS suggests that the results of such 
studies will depend strongly on the choice of discretization, particularly in the discretization of the nonlinear 
term. Therefore, the choice of discretization for VNLS merits closer analysis. 

The system 

d I 2 

i— q«= Toiqn-i -2q„ + q n+ i)+ ||q n || (q„-i + q„-i) (5) 
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where q„ is an V-component vector, manifests properties if an integrable system (cf. Section ^|). Moreover, 
this system (^|) is extremely useful as a numerical scheme for VNLS (Q). Therefore, the underlying chaotic 
nature of some periodic solutions of (^) and (||) can be examined effectively. Because (||) is a natural 
generalization of IDNLS (^) that has some important symmetries of of VNLS , we refer to this system as 
the symmetric discretization of VNLS or, simply the symmetric system. 

To date, no compatible linear operator scattering pair for the symmetric system ^ has been derived. 
Without such a pair (|) cannot be solved via the 1ST. However, under the reduction q„ = e 4 "*v„, where v„ 
is independent of t, (p) reduces to an V-dimensional difference equation which is known to be integrable (cf. 
p2| ). Also, we note that the reduction q„ = cq n , where c is a constant complex unit vector (||c|| ~ 1) and 
q n a scalar, (||) reduces to the integrable discretization of NLS @. 

The symmetric discretization (||) has a class of analytical multi-soliton solutions || ^o). These exact 
solutions of the symmetric system only reduce to a subset of those associated with (|j) . In order to determine 
whether general multi-soliton solutions exist, we examined, by numerical simulation, the interactions of 
solitary waves associated with the symmetric system which lie outside this class of known analytical solutions. 
The numerical evidence, discussed later in this paper, indicates that, in the symmetric system, the solitary 
waves interact elastically and are therefore true solitons. This finding strongly suggests that the symmetric 
system (||) is indeed integrable. 

In the future, we shall identify of regions of phase space where effective chaos could be expected for the 
symmetric system. In these regions of phase space, the dynamics of the symmetric discretization can be 
compared, by numerical simulation, to the dynamics of different discretizations of VNLS- e.g. replacing the 
nonlinear term in (0) by ||q„|| 2 q,i in analogy with (H). 



2 Effective Chaos in NLS 

Let us return briefly to the scalar equation (|l|). In both discretizations (|^,|3|), truncation error is 0{h 2 ). 
For the non-integrable discretization (Q) the truncation error is a perturbation from integrability while the 
integrable discretization (j^) is, as its name implies, integrable for finite h. For NLS with periodic initial 
data, some initial conditions evolve irregularly in space and time when numerically integrated via DDNLS 
at moderate discretizations. At high levels of discretization this irregularity disappears though often 
very slowly. This "chaotic" response is due to the perturbation from integrability in the truncation error in 
DDNLS (||) , the non-integrable discretization. Simulations with IDNLS ([|) and the same initial data do not 
show irregularity even at moderate levels of discretization EL |l^| . 

Irregular and chaotic dynamics can also be caused by perturbations due to round-off error in the floating 
pont arithmetic of computer simulation. This is demonstrated by integrating the initial data with the 
integrable discretization ([|) via the mathematically equivalent discretization 



i — q n = To {Qu-1 - 2 <7n + q n +l) + |gw| 2 <Zn+l + |<Zn| 2 <?n-l- (6) 
At h A 

The only difference between the discretizations (^) and (^|) is the distribution of the multiplication in the 
nonlinear term. Nonetheless, in numerical simulations, the two solutions, q n (t) and q n (t), with the same 
initial data (q n (0) — <?n(0)) can dramatically drift apart until their difference is comparable to the amplitudes 
of the solution || fjj. 

These instances of chaotic behavior in NLS are associated with proximity of the initial data to homoclinic 
manifolds of NLS. The location and structure of these homoclinic manifolds can be determined for NLS by 
the spectral theory of the associated linear scattering problem. When initial data are chosen well away from 
the homoclinic manifolds, the chaotic responses described above are observed only at coarser discretizations 
or not at all In low dimensional systems, chaotic dynamics in regions of phase space are associated 

with crossings of the homoclinic manifold. However, by analyzing the spectral data of the associated linear 
problem, one can demonstrate that for DDNLS, chaotic dynamics generically occur without the solution 
crossing the homoclinic manifold S . 
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3 Vector NLS 



Solutions of VNLS (|J) are more complex than those of the scalar case ([!]) because there are more degrees 
of freedom. In particular, the dynamics of soliton interactions depends on the vector nature of the system. 
A clear definition of solitons in the vector system is a prerequisite for analyzing soliton solutions in the 
discretizations of VNLS such as the symmetric system (||). Hence, here, we give a brief account of some key 
properties of VNLS and their affect on the soliton solutions. 
Under the reduction 

q = cq (7) 

where c is a constant, iV-component vector such that ||c|| 2 = 1 and q is a scalar function of x and t, vector 
NLS (|) reduces to NLS (@). This is a manifestation of the fact that the vector system (J4J) is a generalization 
of NLS. When the vector- valued independent variable q is restricted to a one-dimensional subspace (which is 
equivalent to a scalar), NLS is recovered. As a consequence, any solution of scalar NLS has a corresponding 
family of solutions of VNLS. We call a solution of the form (^) a reduction solution. We refer to c as the 
polarization of the reduction solution. 

The solitons of VNLS can be found by 1ST for the boundary condition of rapid decay on the whole 
line The one-soliton solutions of VNLS are the reduction solutions (i.e. in the form of eq. [?]) in 

which the scalar function q is a one-soliton solution of NLS and c can be any unit vector. In the backward 
(t — » — oo ) and forward (t — > +oo) long-time limits, a generic rapidly-decaying whole-line solution of VNLS 
asymptotically approaches a linear superposition of solitons: 



\ + 



■j 



2^: 



as t — > ±oo 



where = cfqf, is a complex unit vector and q^ is a one-soliton solution of NLS such that has the 
same amplitude and speed as q^ [jT6| . Therefore, in the long-time limits, we can define a separate polarization 
vector for each sc 
that, in general, 
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vector for each soliton, namely c^. Comparison of the forward (+) and backward (— ) long-time limits shows 



The polarization of each individual soliton shifts due interaction with the other solitons. However, 

||c+|| 2 = ||cT|| 2 = l. (8) 

In addition to the shift in the polarization, the locations of individual soliton peaks are shifted by the 
interaction with other solitons (this is the usual phase shift seen in scalar soliton equations). A closed 
formula for these phase shifts (the change in polarization and in the location of the peak) can be derived via 
the 1ST H| . 

We refer to the squared absolute value of a component of the polarization vector as the intensity of that 
component of the polarization. The intensity of the ^-th component of the j-th soliton is given by \c*p + \ 2 

(\cj | 2 ) where c"p + (cp ) is the £-th component of cj~ (cj), the polarization vector of j-th soliton in the 
forward (backward) long-time limit. The relation (||) implies that, for each soliton, the sum of the intensities 
is constant, namely equal to one. However, the distribution of the intensity among the components of each 
soliton will not, in general, be equal in the forward and backwards long-time limits due to interaction with 
other solitons: that is, subject to the constraint of eq. (||), in general, 

|cf -|V kf + | 2 - (9) 



4 



This shift in the distribution of intensity between the backward and the forward long-time limits (|]) is a 
distinctive feature of soliton interactions in the vector system (^). There is no corresponding phenomenon 
in the scalar equation ([!]). 

Manakov showed by 1ST that the following special case holds for Af-soliton solutions: if, for any 
j,k = l ... M, either 



(10a) 



or 



|ct. c -| = 0, (10b) 
where • denotes the inner product, then it turns out that 

|C7. C +| = 1 (11) 

and, moreover, 

| C f-|2 = |c W +| 2 (12) 

for all £ = 1, 2, j = 1, . . . , M (we assume the case of two components). These soliton solutions constitute a 
special class of vector solitons where the distribution of intensity of an individual soliton is not shifted (i.e. 
the relation ([12]) holds) even though the solitons pass through each other. 

When considered as a model for electromagnetic wave propagation in optical fibers, the components of 
VNLS (Q) play the role of a basis for the polarization vector of the transverse electric field. Because this 
choice of basis is arbitrary, the vector system (0) ought to be, and is, invariant under a change of basis. 
Mathematically, a change of basis is obtained by the independent- variable transformation 



Uq 



(13) 



where U is a unitary matrix. The new independent variable, q, satisfies VNLS if, and only if, q does. We 
note that, in particular, the distribution of intensity for a soliton depends on the choice of basis. However, 
the conditions (10a- 10b) and the consequence (11) are independent of the choice of basis- i. e. dot products 
are preserved under unitary transformations. 



4 Discretization of Vector NLS 

As we noted previously, the symmetric system (|5|) is the natural vector generalization of IDNLS (||). However, 
the most natural vector generalization of the linear scattering pair for IDNLS yields the following nonlinear 
equation as a result of compatibility: 



*^ q ™ = 7^ ( q ™~ 1 ~~ 2q " + q " +1 ) ~~ ( r « £ l»- 1 ) C l» ~ ( r n c ln) c ln+l 



(14a) 
(14b) 



where q„ 



is an TV-component vector as is r„ and the superscript T denotes the transpose. To get VNLS 
q* (* denotes the complex conjugate) after taking the continuum limit. The 



from (14aT4b), we take r 

discrete system (14a- 14b| ) does not admit the symmetry r n = — q* for finite h. Hence, we refer to ( 14a - 141 j ) 
as the asymmetric discretization of VNLS, or more briefly as the asymmetric system. We emphasize that 
(14a 44b) is not equivalent to the symmetric system (||) 
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The asymmetric system has a class of soliton solutions which corresponds to all the soliton solutions of 
the continuous version of VNLS (^) ||. However, in these solutions, |r„ — (— q*)| = 0(h) so the asymmetry 
cannot not be ignored In general, for initial data such that r„ = — q*, the solution will evolve so that 
r n =/= — q* at later times. The unavoidable asymmetry of the asymmetric system (14a-14b) makes it less 
desirable as a discretization of VNLS. Hence, we turn our attention to the symmetric discretization. 

We now discuss some important properties of the symmetric discretization. The symmetric system (^) 
reduces to IDNLS (|^) under the reduction 



q« = cq n (15) 

where c is a constant vector such that ||c|| = 1. This is the discrete analog of eq. (0). As for VNLS, we 
refer to the vector c in ( |l5|) as the polarization of the reduction solution q„. Also, the symmetric system is 
invariant under a discrete version of the independent-variable transformation (^). Hence, as for the PDE 
, the components of q„ can be considered to be a basis of the polarization and the discrete form of 
the unitary transformation ( |l3| ) as a change of coordinates. Perhaps the most important property of the 
symmetric system is the existence of soliton solutions. These are discussed in the next section. 



5 Soliton Solutions of the Symmetric System 

In order to demonstrate the existence of solitons in the symmetric system, we must first identify the appro- 
priate vector solitary waves. Given these solitary wave solutions, we then investigate their interaction. To 
construct a vector solitary- wave solution of the symmetric system (^|) multiply the scalar one-soliton solution 
of IDNLS (0) by an arbitrary unit vector. The resulting vector solution q„ is the reduction solution- i.e. of 
the form (p5|) -where the scalar function q n is given by: 



q n (t) = Ae-^-^-'t'hecHahn - vt - 0) = e^q n {t) (16) 

with 



sinh(ah) 2(1 — cosh(ah) cos(bh)) sinh(a/i) s'm(bh) 

A= ~h~' " = h* ' V = 2 h? 

where a,b,9,(j) are arbitrary parameters. For this vector solitary wave, a and b determine the amplitude, A 
and speed, v. The parameter 9 is the position of the peak of | |q| | at t = and 4> is an overall complex phase. 

The simplest example of the interaction of these vector solitary waves is constructed from the scalar 
M-soliton solution of IDNLS. Let qj. n be the scalar M-soliton solution of IDNLS with parameters aj, bj 
and consider the vector solution q n = cq n where, as before, c is a unit vector. Then in the long-time limits 
(t -> ±oo), 



M 

q±~cX>^«t»' ( 17 ) 

3=1 

where qp n is a solitary- wave solution of the form given above in eq. ( |l6[ ) with a — aj, b = bj, 9 = 8^. The 
parameters (jy^ are complex phase shifts of the scalar soliton interaction. For such a solution, each solitary 
wave envelope can be thought of as having an individual polarization: 



= ce^t . 

In this example, we see that each of the vector solitary waves has the same amplitude and speed in the 
forward long-time limit as it had in the backward long-time limit (i.e. a soliton- type interaction). This 
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suggests that the individual vector solitary waves constructed from ( |l6| ) are, in fact, the solitons of the 
symmetric system (|^). 

However, in this simple example (|l^), since c is the same for each solitary wave, the polarizations of the 
vector solitary waves must satisfy 

l c 7 • c aTI = i 

for all j,k = 1 . . . M. There is no such constraint for the solitons in the vector PDE (^). Therefore, we 
investigate the more general vector solitary-wave interactions (i.e. \cj • | < 1) in the discrete symmetric 
system by direct and numerical methods. 

5.1 Soliton solutions by Hirota's Method 

To find soliton solutions of the symmetric system (||) by Hirota's method we write it in the bilinear form: 
by taking 

gn 

q« = -r- 

Jn 

where g„ is a vector and /„ is a scalar, we find 

ih 2 D t f„ ■ g„ = /n-lgn+l/n-l - 2/ n g„ + /n+lgn-l (18a) 
fn+lfn-1- f^h 2 ^ 2 (18b) 

A two-soliton of the symmetric system is given by the following solution of the bilinear equations (|l8af- 



18b): 



/„ = 1 + e ^- n+7 i*^ + e^+vi,i + \B 1 \ 2 e lll ^ +r ''.n+m,n+vln (ig a ) 

" } = \ ( ehPl ~ e ~ HP ' 1 ) eVl '" i 1 + B ^ m ' n+nln ) ( 19b ) 
W = I (e hp2 - e-' lp '^ e" 2 - (l + B 2 e^»+^.») (19c) 



where is the k-th component of g„ and 



Vi,n = Pjnh + ± (2 - e hp * - e- hp t ) t. (20) 



The coefficients Bj are 



( e h Pi — e h P2^ e hpi _j_ e hP2\ (e hpi — e hp2 )(e hp * 1 + e hp2 ) 



! e h{pi+p 2 ) -|- l)( e M>i+Pa) — 1) ' ( e h(pi+P2) -|_ l)( e ^bi+w) — 1) 

and the complex constants pj — cij — ibj (where a,j > and j = 1, 2) determine the amplitudes and envelope 
speeds of the solitons. 

In order to see that the above solution ( 19a - 19c| ) indeed gives a two-soliton solution, consider the long- 



time limits (i.e. t — > ±oo). In the forward (t — » +oo) and backward (t — > — oo) long-time limits, this solution 



(19a-19c) asymptotically approaches the form 



~ qr^n + q^n 
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where q~ n — c f<lj n an d j = 1,2. The vectors cj are unit vectors and the scalars q^ n are solitary waves of 
the form ( ^6| ) with a = aj, b = bj and 6 = 9^ where j = 1,2. If, without loss of generality, we assume that 
cij and 6j are such that sinh(oi/i) sin(6i/i) > sinh(a2/i) sinfah) then 



and 



where 



?r = 0, 0+ = -log|£i|, 02 = -logical, ^ 2 + = 



cf = (X,0), 



if = + arg Si , 



hbi + arg _B 2 , 



hb 2 



The two-soliton solution (19a- 19c) does not encompass the most general case because 



(21a) 



and 



= 0. 



(21b) 



This solution can be multiplied by a unitary matrix (the discrete analog of the transformation in eq. |l3| ) to 
obtain a two-soliton solution with any pair of vectors such th at ([21a|) holds (under such a transformation, 
(21h) will still hold). Therefore two-soliton solutions given by fll94|l9c] ) are a restricted class of two-soliton 
solutions (where eq. (21a- 21bj ) constitutes the restriction). 

For more than two solitons, Hirota's method can still be applied. The M-soliton solutions of the sym- 
metric system take on a considerably more complicated form; they can be represented as ratios of Pfafhans 
lEOfl. Sill, for these M-soliton solutions, cither 



= 1 



(22a) 



< >r 



(22b) 



and |cj • ct| = 1 for j,k = 1 . . . M @ |. 



Even though analytical formulae exist for multi-soliton solutions of the symmetric system (with any 
number, M, of solitons), these solutions are constrained to satisfy (22a- 22b| ). Such solutions are more 
general, but are like th e special class of soliton solutions of the PDE discussed at the end of Section |[ They 
satisfy the restriction (21a -21b). 



5.2 General Soliton Interactions by Numerical Simulation 

Recall that, in the PDE <Q), the polarizations of the individual solitons in the multi-soliton solutions (derived 
by 1ST) are not constrained- i.e there are vector multi-soliton solutions of the PDE for which the polarizations 
satisfy < \cj -c^ \ < 1 for some j, k. For the discrete symmetric system (|H|), we investigated the general 
solitary-wave interaction- i.e < |cj ■ cj^| < 1 by numerical simulation because, to date, there are no 
analytical formulae covering this general case. Moreover, it is important to determine whether eq. @ is 
integrable. The existence of general soliton solutions would strongly support the possibility that (ra) is indeed 
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integrable. In addition, as we shall show, the measurement of soliton parameters in discrete systems requires 
a novel approach. 

In these simulations, the initial condition was composed of vector solitary waves q~ n of the form qj n — 
cj qj n where the scalar qj n is of the form (^|). In the initial condition, these vector solitary waves were well- 
separated. Then, the symmetric system was integrated in time by an adaptive Runge-Kutta-Merson routine 
(from the NAG library) until the peaks were again well-separated (see Figure [j]) . The separation of peaks in 
the initial and final conditions makes these finite-time conditions comparable to, respectively, the backwards 
and forwards long-time limits. In the example Figure [l], and in other simulations, the solitary waves appear 
to interact without any any radiation or loss. We confirmed this visual result with striking quantitative 
measurements: we measured the difference between the solitary waves resulting from the numerical simulation 
and the exact solitary wave form, 



%> = 4Un ( 23 ) 

where q- n is of the form (^6|). 

In order to compare the final-time data of the numerical simulation with exact solitary-wave solutions it 
is necessary to estimate the polarization, c" 1 " , and the shift in the envelope, 6j for the j-th solitary wave at 
the final time (the values of cj, and 8 J are fixed by the choice of initial data). If the j-th solitary wave is 
actually of the form ( p3| ) at the final time, then the £-th component of the vector satisfies 



„W+ _ i{bjhn~u>jtf) _ 



1j,n 



kS + l 2 + lCTl 



(2)H 



(24) 



for all n where ||qj~ n || is not negligible (we denote this set of points n containing the j-th solitary wave by 
fL). In practice, we compute the estimate c\ by 



avg„ e o 3 



Jl),/|2 
q-a I 



(2),/| 



(1) f f2l f 

where £ = 1,2 and (q„ , qn ) is the numerical data at the final time, t = tf. Similarly, if the absolute 
value of j-th solitary wave is a sech envelope as in (pi) , then for all n G , 



,(1)+|2 



9^ = ahn — Vjtf — sech 1 
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(2)+|2^ 



A, 



(25) 



where Aj 



sinh(a jh) 



Again, we construct the estimate 8f by taking the average of the right-hand side of 



h . ^ „.„, „„„ „„„„„„„„ „^ 

( |25| ) over n G ilj where we substitute {qf n + ,qf n + ) with (qiP J ,q { n )J ), the numerical data at the final time. 
Note that we assume that aj and bj do not change in the evolution and hence we do not estimate them in 
the forward long-time limit. 

The difference between the numerically-derived solitary wave and a solitary wave-form of the type 
(p3|) is measured by 



max 



A, 



(26) 



where q 1 



(qn^'^ : qn^'^) for n e flj and q^ n is the estimated soliton wave-form 



-A j e- i to ka - u ' t risech (djhn - v 3 t f - 8+^j 
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In order to study the symmetric case in a parameter regime that is "far" from the continuum limit, we ran 
the simulations with amplitude parameters, Oj, and grid sizes h such that the soliton width was comparable 
to the grid size (equivalently, a>jh w 1). For the parameters bj we considered two regimes: (i) bj such that 
the spatial frequency of the complex carrier-wave e ibjnh was comparable to the grid size (bjh ~ 1) and (ii) 
bj such that the spatial frequency was large compared to the grid size (bjh .1). We ran simulations with 
several values of cij for each bj. In particular, when bjh « .1 the speed of the solitary-wave envelopes was 
slow, the solitary waves interacted over a long time. 

The time-integration routine we used (routine D02BAF from the NAG library) is adaptive and computes 
until an internally-defined error tolerance reaches a user-specified value. In the simulations, when the error 
tolerance in the adaptive integration scheme was decreased, the difference between the data at the final time 
and a soliton waveform (26) decreased proportionally. Table [j] shows data from an example experiment. In 



this experiment, the initial data is generic in that < |c^ • c^~| = .6 < 1 and, to date, there is no known 
analytical formula for the solution. The error, A,-, is proportional to the error tolerance of the integration. 
This indicates that the numerically computed solution converges to an elastic soliton interaction as the time 
integration is computed more accurately. We conclude from numerous such experiments that, in fact, the 
vector solitary waves of the form qj n = c~qj n (where the scalar q~ n is of the form of eq. [if]) interact 
elastically- i. e. the vector solitary waves behave as true solitons -in the symmetric system (||). 

For a generic interaction of two solitary waves of the symmetric system, < |c^ ■ c^~| < 1. In the PDE 
(|J) , such generic vector soliton interactions result in the shift of the polarizations of the individual vector 
solitons- i.e. \cj ■Cj] < 1, j = 1, 2. We observe the same distinctive vector behavior in numerical simulations 
of the discrete symmetric system. Table shows the results of a number of typical experiments. In a two 

. — (2) — 

soliton-interaction we can always pick a basis for the polarization such that c l = (1,0) and argCj = 0. 

— . . (2) — 

Because c 2 is a unit vector with two components, the absolute value of c 2 satisfies: 



l4 2) ~l = V 1 

Hence, as in Table ^ in order to consider a general two-soliton interaction, we need only vary ■ In our 
simulations, we see the polarizations of the individual solitary waves shift due to interaction while the solitary 
waves retain their shape. In the example interactions tabulated in Table ^, there are initial conditions that 
result in large polarization shifts, that is the inner products of the forward and backward polarization of an 
individual vector solitary wave is small (e.g. |c^ • c^| = .158). This is the intrinsically vector-type soliton 
interaction of the solitary waves. Therefore, we conclude from the numerical simulations that the symmetric 
system displays general vector-soliton behavior analogous to VNLS, the PDE continuum limit. 
We also considered the following generalization of the two-component symmetric system 

= - 2^ + «&) + (l^l 2 + B\ q ^) (9& + ffSx) (27a) 

i±£> = ^(A - + &) + I 2 + ki 2) I 2 ) (& + A). (27b) 

In the case \B\ ^ 1, we found very strong inelastic effects. In particular, depending on the initial soliton 
parameters, numerical simulations of soliton interaction, the solitons (i) produced large radiative tails and/or 
(ii) reflected from one another. These numerical results suggest that the general elastic interactions obse rved 
in th e symmetric system (||) are a manifestation of the integrability of that system. We mention that ( |27a| - 
27b ) is a an example of the interesting generalization of where the nonlinear term for the j-th component, 



(Ef=i \Qn ] \ 2 ) + is replaced by (£f =1 | 2 ) (g^ + q^ +1 ). We expect that this model 

will exhibit a range of different inelatic effects between solitary waves such as those observed in (27a- 27h| ). 

The mechanism of the elastic soliton interactions observed for the symmetric system remains to be 
completely explained analytically. More generally, a proof that the symmetric system can be completely 
integrated remains as an important open problem. 



10 



Acknowledgments 



This work was supported in part by NSF grants DMS-9404265 and DMS-9703850 as well as U. S. Air Force 
AASERT grant DOD AF-F49620-93-I-0574 and U. S. Navy AASERT grant DOD N00014-94-1-0915. 

References 

M. J. Ablowitz, S. Chakravarty, and B. M. Hcrbst. Intcgrability, computation and applications. Acta 
Applicandae Mathematicae, 39:5-37, 1995. 

M. J. Ablowitz and P. A. Clarkson. Solitons Nonlinear Evolution Equations and Inverse Scattering. 
Number 149 in London Mathematical Society Lecture Note Series. Cambridge University Press, 1991. 

M. J. Ablowitz, B. M. Hcrbst, and C. M. Schobcr. Computational chaos in the nonlinear Schrodingcr 
equation without homoclinic crossings. Physica A, 228:212-235, 1996. 

M. J. Ablowitz and J. F. Ladik. A nonlinear difference scheme and inverse scattering. Studies in Applied 
Mathematics, 55:213-229, 1976. 

M. J. Ablowitz, Y. Ohta, and A. D. Trubatch. On discretizations of the vector nonlinear Schrodinger 
equation. APPM Preprint 349, University of Colorado-Boulder, March 1998. 

M. J. Ablowitz and C. M. Schober. Effective chaos in the nonlinear schrodinger equation. Contemporary 
Mathematics, 172:253-268, 1994. 

M. J. Ablowitz, C. M. Schober, and B. M. Herbst. Numerical chaos, roundoff errors and homoclinic 
manifolds. Physical Review Letters, 71(17):2683-2686, October 1993. 

M. J. Ablowitz and H. Segur. Solitons and the Inverse Scattering Transfom. Number 4 in SIAM Studies 
in Applied Mathematics. SIAM, 1981. 

D. Cai, A. R. Bishop, N. Gr0nbech- Jensen, and M. Salerno. Electric field-induced nonlinear bloch 
oscillations and dynamical localization. Physical Review Letters, 74:1186, 1995. 

J. C. Eilbeck, P. S. Lombdahl, and A. C. Scott. The discrete self-trapping equation. Phyisica, 16 
D:318-338, 1985. 

S. G. Evangclidcs, L. F. Mollcnaucr, J. P. Gordon, and N. S. Bergano. Polarization multiplexing with 
solitons. Journal of Lightwave Technology, 10(l):28-35, January 1992. 

B. M. Herbst and M. J. Ablowitz. Numerical chaos, symplectic integrators and exponentially small 
splitting distances. Journal of Computational Physics, 105:122-132, 1993. 

V. M. Kenkre and D. K. Campbell. Self-trapping on a dimer: time-dependent solutions of a discrete 
nonlinear Schrodingcr equation. Physical Review B, 34:4959-4961, 1986. 

V. M. Kenkre and G. P. Tsironis. Nonlinear effects in quasilinear neutron scattering: exact line calcu- 
lation for a dimer. Physical Review B, 35:1473-1484, 1987. 

T. I. Lakoba and D. J. Kaup. Perturbation theory for the manakov soliton and its applications to pulse 
propagation in randomly birefringent fibers. Physical Review E, 56(5):6147-6165, November 1997. 

S. V. Manakov. On the theory of two-dimensional stationary self-focusing of electromagnetic waves. 
Soviet Physics JETP, 38(2):248-253, 1974. 

D. W. McLaughlin and E. A. Overman II. Surveys in Applied Mathematics, volume 1, chapter Whiskered 
Tori for Intcgrable Pde's: Chaotic behavior in near Integrable Pde's. Plenum Press, New York, 1995. 

C. R. Menyuk. Nonlinear pulse propagation in birefringent optical fibers. IEEE Journal of Quantum 
Electronics, QE-23(2):174-176, February 1987. 



11 



[19] C. R. Menyuk. Pulse propagation in an elliptically birefringent Kerr medium. IEEE Journal of Quantum 
Electronics, 25(12):2674-2682, December 1989. 



[20] Y. Ohta. Pfaffian solutions for coupled discrete nonlinear schrodingcr equation. Chaos, Solitons and 
Fractals, to appear. Proceedings of Brussels Meeting II: Integrability and Chaos in Discrete Systems 
(Brussels, 2-6 July 1997). 

[21] V. S. Shchesnovich and E. V. Doktorov. Perturbation theory for solitons of the Manakov system. 
Physical Review E, 55(6):7626-7635, June 1997. 

[22] Y. B. Suris. A discrete-time Gamier system. Physics Letters. A, 189(4):281~289, 1994. 



12 



t = -14. 



t = 3.63 



2.98 - 



2.98 - 



1.49- 



1.49- 



. 

t = -8.14 



18 . 



A 



. 

t = 9.52 



18 . 



2.98 - 



2.98 - 



1.49 




1.49- 




. 

t = 15.4 



18 . 



2.98 - 



2.98 - 



1.49- 




1.49 



. 



18 . 




Figure 1: Two-soliton interaction for the symmetric system. The filled boxes are \qn^\ and the open boxes 
arc \qk \- Increasing time is read down column- wise. Soliton 1 is on the left and Soliton 2 is on the right 
at t = —14. They are reversed at t = 15.4. The soliton parameters are: a\ = 1, <Z2 = 2, b\ = .1, 62 = — .1. 
The polarization vectors are: = (1,0), = (.6e iTr/3 , .8) , at t = -14, and cf = (.16e i ' 587r , .98e^- 167r ), 
= (.82e l 277r , .57e~ 4 - 037r ) at t = 15.4. This is a typical two-soliton interaction where < |cj~ • = -6 < 1, 



of I = .16 < 1 and |c 9 



.+ 1 - 



= .95 < 1. 
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Error: log 10 Aj 


tol 


soliton 1 


soliton 2 


ict" 


-4.11 


-5.97 


10- 7 


-5.33 


-6.99 


10- 8 


-6.58 


-8.00 


10- 9 


-7.84 


-9.01 


10 -io 


-9.12 


-10.0 



Table 1: Difference between solution at the final time and a soliton wave- form for a two solitary- wave 
interaction. Soliton parameters: a\ = 2, b\ = 0.1, = (1, 0), 0,2 = 1, b% = —0.1, = {-^= -^=). "tol" 
is the error tolerance allowed by the NAG routine D02BAF. The errors in the table are the log 10 of the 
error given by (|26| ) . The error decreases proportionally with the error tolerance used in the simulation. This 
indicates that the errors are due to the time integration and that, therefore, the solitary waves interact as 
solitons. 



C 2 


|cr-cri 


\ C 2 - C Z\ 


C 2 


|cr-cj-| 


l C 2 ' C 2~ 


(0.2e i0 , .98e i0 ) 


0.843 


0.985 


(0.2e l7r / 3 , .98e l °) 


0.843 


0.985 


(0.4e i0 , .92e i0 ) 


0.465 


0.958 


(0.4e l7r / 3 , .92e l °) 


0.465 


0.958 


(0.6e l0 ,0.8e io ) 


0.158 


0.947 


(0.6e M / 3 ,0.8e i0 ) 


0.158 


0.947 


(0.8e*°,0.6e i0 ) 


0.544 


0.962 


(0.8e l7r / 3 ,0.6e l °) 


0.544 


0.962 


(0.2e"/ 2 ,.98e io ) 
(0.4e"/ 2 ,.92e io ) 
(0.6e i7r / 2 ,0.8e l °) 
(0.8e iTr / 2 ,0.6e io ) 


0.843 
0.465 
0.158 
0.544 


0.985 
0.958 
0.947 
0.962 


(0.2e in , .98e l °) 
(QAe™, .92e l °) 
(0.6e l7r ,0.8e l °) 
(0.8e l7r ,0.6e") 


0.843 
0.465 
0.158 
0.544 


0.985 
0.958 
0.947 
0.962 



Table 2: Shift in polarization in the two-soliton interaction of the symmetric system. The soliton parameters 
are a\ = 1, 02 = 2, b\ = .1, 62 = — .1. The polarizations before interaction are = (1,0) and as 
given in the table. The polarizations after interaction are cf and cj. The values \cj • cj~| < 1 indicate 
that the polarization vectors are shifted by the soliton interaction. These results were obtained by numerical 
simulation. 
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